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1 INTRODUCTION 



The observational study of the early stages of galaxy for- 
mation is starting a golden age. Among of the best tar- 
get populations are galaxies with strong emission in the 
Lyman-Q line, known as Lyman Alpha Emitters (LAEs) 
(|Partridge fc PeebleiJ[l967l ). Large observational samples of 
these galaxies at redshifts 3 < z < 7 are already avail- 



ABSTRACT 

Using CLARA (Code for Lyman Alpha Radiation Analysis) we constrain the escape 
fraction of Lyman-a radiation in galaxies in the redshift range 5 < z < 7, based on 
the MareNostrum High-z Universe, a SPH cosmological simulation with more than 2 
billion particles. We approximate Lyman-a Emitters (LAEs) as dusty gaseous slabs 
with Lyman-a radiation sources homogeneously mixed in the gas. Escape fractions 
for such a configuration and for different gas and dust contents are calculated using 
our newly developed radiative transfer code CLARA. The results are applied to the 
MareNostrum High-z Universe numerical galaxies. The model shows a weak redshift 
evolution and good agreement with estimations of the escape fraction as a function 
of reddening from observations at z 2.2 and z ^ 3. We extend the slab model by 
including additional dust in a clumpy component in order to reproduce the UV con- 
tinuum luminosity function and UV colours at redshifts z > 5. The LAE Luminosity 
Function (LF) based on the extended clumpy model reproduces broadly the bright 
end of the LF derived from observations at z ~ 5 and z ~ 6. At z ~ 7 our model 
over-predicts the LF by roughly a factor of four, presumably because the effects of the 
neutral intergalactic medium are not taken into account. The remaining tension be- 
tween the observed and simulated faint end of the LF, both in the UV-continuum and 
Lyman-a at redshifts z ~ 5 and z 6 points towards an overabundance of simulated 
LAEs hosted in haloes of masses 1.0 x IO^^/i-^Mq < Mh < 4.0 x IQ^^H-^Mq. Given 
the difficulties in explaining the observed overabundance by dust absorption, a prob- 
able origin of the mismatch are the high star formation rates in the simulated haloes 
around the quoted mass range. A more efficient supernova feedback should be able to 
regulate the star formation process in the shallow potential wells of these haloes. 

The importance of LAEs is not only limited to galaxy 
evolution. The detailed measurement of their clustering 
properties, in pa rticular the Baryonic A coustic Oscillations 
(BAO) feature (|Eisenstein et al.1l2005l ). is expected to be 
detected at high redshifts, potentially providing useful con- 
straints on the evolution of dark energy. A key point in this 
analysis is understan ding the bias of LAE s as tracers of the 
large scale structure (|Wagner et al.|[2008l ). 
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J iected to be ga thered in ongoing and future observations 
Hill et al.lliooa '). 



The study of the epoch of reionisation has also greatly 
benefited from the study of LAEs, not only because the fea- 
tures of the emission line make its observational detection 
unambiguous, but also because the Lyman-Q photons are 
sensitive to the distribution of neutral hydrogen, and the 
changes in the line are able to constrain the ionisation state 
of the IGM. It is thus of crucial importance to properly 
model the propagation of Lyman-Q photons through arbi- 
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trary gas distributions which might c ontain some dust c aus- 
ing the absorption of these photons (|Zheng et aLll2010l ). 

Because of the resonant nature of the hue, Lyman-Q 
photons perform a random walk in space and frequency be- 
fore esc aping the neutra l ISM/IGM and reaching the ob- 
server iHarringtonl 1 19731 ). As a result, the sensitivity of 
a Lyman-a photon to dust absorption is enhanced. Small 
quantities of dust, depending on the amount and dynamical 
state of the neutral gas, can significantly diminish the in- 
tensity of the Lyman-a line. The detailed shape of the line 
profile also depe nds on the dy namical state of the gas and 
its dust content (|Neufeldlll9"9ll '). 

The observational estimation of the fraction of Lyman- 
a photons escaping the ISM (hereafter escape fraction) is 
challenging. It usually requires another probe (UV contin- 
uum or a non-resonant recombination line such as Hq) and 
some estimation of the continuum dust extinction. Recent 
constraints at 2: ~ 2.2 are based on blind surveys of Lyman- 
a and Ha. As the line ratio between Ha and Lyman-a is 
constant and given by atomic physics, the measurements of 
the Ha line intensity corrected by extinction, allow for the 
estim ation of the intrinsic Lyman-a emission es et al.l 
|2010D . 

Concerning the emission process, there is general agree- 
ment that the bulk of Lyman-a luminosity in high redshift 
LAEs is triggered by star formation processes. The fraction 
of the emission coming from coUisional excitation of the gas 
is small and there is evidence th at AGN activities do not 
power the Lyman-a line emission jWang et al.ll2004l ). 

Due to all the efforts to model LAEs in a cosmological 
context, it has been recognised that the predicted abundance 
of LAEs (when considering the intrinsic emission) overes- 
timates by orders of magnitude the obs erved abundance 
as quan tified by the l uminosit y function ([Le Delliou et al.l 
[2005; K obavashi et al.|[2007, : . Zheng et al.ll201(]h . Asthenum- 
ber density of host dark matter haloes is tightly constrained 
by the allowed range of cosmological parameters in the 
ACDM cosmology, the most reasonable assumption is that 
the difference between observed and predicted abundances 
is due to the absorption of Lyman-a photons in the ISM. 
These facts motivate the crucial importance of a theoretical 
determination of the escape fraction of Lyman-a emitting 
galaxies. 

The theoretical models of LAEs populations based 
on numerical simulations fall into two families: (i) semi- 
analytic models an d (ii) hydrodynamical models. On the 
semi-analytic side, iLe Delliou et al.l (|2005l ) assume a con- 
stant escape fraction for all galaxi es, with values c l ose to 
fesc ~ 0.02. In a more refined model iKobavashi et al.l (|2007l ) 
allow a variable escape fraction motivated by a wind feed- 
back model. In the realm of hyd rodynamical models the re- 
cent work of IZheng et al.l (|201G| ) addresses the effect of the 
IGM with full radiative transfer of the Lyman-a line. They 
find that as a consequence of pure photon diffusion a LAE 
can show a low surface brightness, and might be missed in 
a survey. This effect introduces an effective escape fraction 
that is not the result of dust absorption. 

Hydrodynamical simulations of single galaxies 
iLaursen et al.l |2009| ) and simp le gas/dust configura - 
tions have been explored as well (jVerhamme et aLll2006l ). 
The problem with simulated individual galaxies is that the 
small sample available so far does not allow for the infer- 



ence of useful scalings or valid statistics in a cosmological 
context. The limitations of studying simplified configu- 
rations is that, even though they allow for a wide range 
of models to be simulated, the parameters, such as dust 
abundance or star formation rates, are not constrained by 
any other assumption, making it impossible to infer possible 
scalings with galactic properties already constrained by 
obs ervations. 

iDaval et al.1 (|201G| ') have used a low resolution cosmolog- 
ical SPH simulation to fix the mass contents and star forma- 
tion rates of the galaxies. Unfortunately, they only consider 
the radiative transfer effects of the Lyman-a as a free pa- 
rameter in their model, and do not attempt to bound the 
escape fraction from physical considerations consistent with 
the resonance nature of the Lyman-a line. 

In this paper, we address the problem of deriving statis- 
tics on the escape fraction of high-z LAEs between redshifts 
5 < z < 7 due to the effects of the dusty ISM and resonance 
scattering of the Lyman-a line within an explicit cosmo- 
logical context. We base our physical analysis of the escape 
fraction for a given dust and gas abundance on the results of 
a new state-of-the-art Monte Carlo radiative transfer code 
called CLARA (Code for Lyman Alpha Radiation Analysis). 
The astrophysical application of these results relies on the 
analysis of a SPH galaxy formation simulation with 2 bil- 
lion particles, the MareNostrum High-z Galaxy Formation 
Simulation. 

We follow the approach of obtaining the escape fraction 
for a single family of models (homogeneous slab with differ- 
ent optical depths of gas and dust), and applying them to the 
galaxies in the simulation. The dust content has been cal- 
culated in the simulation by matching the behaviour of the 
UV continuum (luminosities and colours) with the obser ved 
estimates at these redshifts (jPorero- Romero et al.|[201ol ). 

Our approach is dictated by two main technical con- 
straints: 1) it is still not feasible to run the radiative trans- 
fer code on several thousands of individual galaxies in the 
simulation box. 2) the mass and spatial resolution in our 
simulation is not high enough for the radiative transfer cal- 
culation to converge on the escape fraction using the gas 
distribution directly from t he SPH simula t ion, a ccording to 
the convergence studies bv lLaursen et all (|2009l ). 

Our theoretical approximation for the gas and source 
distributions is based on the premise of consistency with 
the model we used for dust extinction, but also with the ob- 
jective of improving the description of simulated Lyman-a 
emitting galaxies by including two features of the absorp- 
tion of the Lyman-a line in galaxies that are commonly ne- 
glected. Namely, we consider: 

a) the absorption enhancement by the gas content due to 
the resonant nature of the line 

b) a spatial distribution of the Lyman-alpha regions in the 
galaxy where the Lyman-a photons are not forced to have 
statistically the same probability of being absorbed, as it is 
the case for centrally emitted Lyman-a photons in a sphere, 
shell or slab configuration. 

This paper is structured as follows. In Section [2] we 
describe the simulation and the galaxy finding technique. 
In Section 13] we describe our method to calculate the spec- 
tral energy distributions for the galaxies in the sample, as 
well as our simplified dust extinction model. We review 
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the UV continuum p r opert ies of the sample as derived in 
I Forero- Romero et ahl (|201(]| ). Our model for LAEs is de- 
scribed in Section|4l together with its implications on the es- 
cape fraction and its implementation into the MareNostrum 
High-z Galaxy Formation Simulation. We discuss the impli- 
cations of our model in Section[Sl Finally, we summarise our 
conclusions in Section[6l All the details regarding the imple- 
mentation of the Lyman-a Monte Carlo radiative transfer 
code CLARA can be found in the Appendix El 



2 COSMOLOGICAL SIMULATION AND 
GALAXY FINDING 

The MareNostrum High-z Universe simulatioi|3 follows the 
non-linear evolution of structures in baryons (gas and stars) 
and dark matter, starting from 2; = 60 within a cube of 
50/i~^Mpc comoving on a side. 

The cosmol ogical parameters u sed are consistent with 
WMAPl data (|Spergel et all 120031 ') t.e. = 0.3, fib = 
0.045, = 0.7, as = 0.9, Hubble parameter h = 0.7, 
and a spectral index n — 1. The initial density field was 
sampled by 1024'^ dark matter particles with a mass of 
moM = 8.2 X W^H'^Mq and 1024^ SPH gas particles with a 
mass of JTigas ~ 1.4 x 10^/i~^Mq. The gravitational smooth- 
ing scale was set to 2 h~^kpc in comoving coordinates. The 
simulation has been performe d using the TREEPM-I-SPH 
code GADGET-2 (|Springel2005l ). Furth er details on the phys- 
ical se tup of the code can be found in iForero- Romero et al.l 
(|2010l ) 

We identify the objects in the simulations using the 
AMIGA Halo FindeiEl ( AHF) which is described in detail in 
iKnoUmann fc Knebd (|2009h . AHF takes into account the 
thermal energy of gas particles during the calculation of the 
binding energy. The halo consists only of bound particles. 
All objects with more than 1000 particles, dark matter, gas 
and stars combined, are used in our analyses. We assume a 
galaxy is resolved if the object contains 200 or more stellar 
particles, which corresponds to objects with > 400 parti- 
cles of gas. This ensures a proper estimation of the average 
gas column densities and star formation rates in the nu- 
merical galaxies, in agreement with recent resolution studies 
(jTrenti et al.ll2010l ). 

The favoured cosmological parameters estimated from 
the analysis of recent CMB data are different from the ones 
used in the simulation. We have included an additional cor- 
rection in the galaxy abundance from the different number 
density of dark matter haloes in the co smology used in the 
simulation (WMAPl, with erg = 0.90 ISoergel et all [ioosl l 
and the val ues favored in more r ecent works (WMAP5, with 
0-8 = 0.796 iDunklev et al.|[2009h . 



3 SPECTRAL MODELING AND UV 
CONTINUUM 

In this work, we use the same spectral model and follow 
the same extinction model to calculate the Spectral Energy 
Distribution (SED) for each galaxy, as described in Section 3 



of lForero- Romero et al.] (|2010l '). The photometric properties 
of galaxies are calculated em ploying the stellar pop ulation 
synthesis model STARD UST jDev ricndt et al.|[l999l l. using 
the methods described in iHatton ct al. (200^- We adopt a 
Salpeter Initial Mass Function (IMF) with lower and upper 
mass cutoffs of O.IM© and I2OM0. 

The SEDs are built from the AHF catalogues already 
described. Each star particle in the simulation represents 
a burst of stars of a given initial mass and metallicity 
evolved at a given age. We have added all the individ- 
ual spectra of each star pa rticle to build UV magnitudes 
(|Forero- Romero et al.|[201ol ). This allowed us as well to im- 
plement a dust extinction model in two different stellar pop- 
ulations distinguished by age as we describe in the next 
paragraphs. However, we do not use these SEDs to estimate 
the production rate of ionizing photons, which is noisy for 
galaxies resolved with less than 100000 particles for physical 
reasons detailed in Section [5.21 

The dust attenuation model parametrises both the ex- 
tinction in a homogeneous interstellar medium (ISM) and 
in the molecular c louds around young stars, following the 
physical model of I Chario t & Fall (,2000l l. The attenuation 
from dust in the homogeneous ISM assumes a slab geometry, 
while the additional attenuation for young stars is modeled 
using spherical symmetry. 

We first describe the optical depth for the homogeneous 
interstellar medium, denoted by ri'^^^(A). We take the mean 
perpendicular optical depth of a galactic disc at wavelength 
A to be 



(Nh) 



2.1 X 10^1 atoms cm" 



where Ax/Av is the extinction curve from iMathis et al.l 
(|l983h . Zg is the gas metallicity, (Nh) is the mean atomic 
hydrogen column density and rj = {1 -\- z)~°' is a, factor that 
accounts for the evolution of the dust to gas ratio at dif- 
ferent redshifts, with a > from the available co nstraints 
based on simplified theo retical models (Ilnou3l2003h and ob- 
servations around z ~ 3 (|Reddv et al.ll2006l l. The extinction 
curve depends on the gas metallicity Zg and is based on 
an interpolation between the solar neighbourhood and the 
Large and SmaU Magellanic Clouds {r = 1.35 for A < 2000 A 
and r = 1.6 for A > 2000A). 

The mean Hydrogen column density is calculated as 



(Nh) = Xi-_ 



Ma 



r atoms cm 



(2) 
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where Xh = 0.75 is the universal mass fraction of Hydrogen, 
Mg is the mass in gas, r j is the radius of the galaxy and rup is 
the proton mass. The radius, stellar and gas masses for each 
galaxy are taken from the AHF catalogues, where we have 
verified that computing {Nh) from the galaxy catalogues 
yields, on average, similar results than integrating the 3D 
gas distribution of the galaxies using the appropriate SPH 
kernel, provided that the galaxies are sampled with more 
than ~ 200 gas particles. 

In addition to the foreground, homogeneous ISM extinc- 
tion, we also model, in a simple manner, the attenuation of 
young stars that are embedded in their birth clouds (BC). 
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Stars younger than a given age, fc, are subject to an addi- 
tional attenuation with mean perpendicular optical depth 



rr(A)=(i-l)rr"(A), 



(3) 



where is the fraction of the total optical depth for these 
young stars with respect to that is found in the homogeneous 
ISM. 

Without any correction, the simulated LPs have a 
higher normalisation than the observational ones, which 
seems to be a general feature for all ACDM hyd rodynamical 
simulations at high redshift (|Night et al.ll200^ ). The excess 
can be caused by two different effects, both possibly acting 
at the same time: the physics included in the simulation giv- 
ing rise to excessive star formation rates, or the intrinsic UV 
ought to be corrected by dust extinction. In this section we 
review the results from the explanation by a dust correction 
based on the physical model described above. 

The simple approximation of a dust optical depth pro- 
portional to the gas column density leads to a reddening 
which scales with the galaxy luminosity. Massive and lu- 
minous galaxies are more extinguished than less massive 
ones. This is in agreement with simila r numerical results 
l|Night et all l2006l : iFinlator et"all 120111) and ob servational 
constraints at redshift z ~ 3 ( Shaplev et al.|[200ll ) . Applying 
such a correction to the data at redshifts 5 < z < 7 cannot 
explain the faint end of the LP, because there is still a large 
overabundance of simulated galaxies with respect to the ob- 
servations. However, a constant reddening E{B — V)^ 0.2 
(assuming a Calzetti law) on all galaxies uniformly dims the 
LF, t hus providin g a good match at the faint end. A clumpy 
ISM l|lnouell2005l ) would be a plausible physical model to 
explain this effect, which also has proven to be effective at 
high redshift in producing an a lmost constant reddening a s 
a function of galaxy luminosity (|Forero-Romero et al.|[201ol ) . 

The results obtained from matching the observed UV 
LF between redshifts 5 < 2: < 7 to the LFs derived from the 
simulation hint that all stars younger than 25 Myr must ad- 
ditionally be extinguished with parameters = 0.01 for red- 
shifts z ~ 5, 6 and ^, = 0.03 for redshift z ~ 7, and a = 1.5 
(|Forero-Romero et al.l I2OIOI ). This means that the young 
stars have an additional extinction — 1, ie. 30 — 100, 
times larger than the extinction associated to the homo- 
geneous ISM. Observational constraints at 2 = locate 
around 1/3 with a wide range of scatter between 0.1 and 0.6 
(|Kong et al.ll2003 ). One possible interpretation for the evo- 
lution of the ^ parameter is that high redshift galaxies have 
a less dust enriched homogeneous ISM than those at low red- 
shift, making the relative contribution of extinction around 
their young stars higher. Within this interpretation, a factor 
of ~ 30 increase in the value of /i (with values of a = 1.5) 
would require that the dust to gas ratio in the homogeneous 
ISM increases between 2 ~ 7 and 2 = by at least a factor 
of ~ 30 X (1 -I- z)" ~ 400, which is feasible under conserva- 
tive theoretical considerations of what the dust to gas ratio 
evolution should be. For instance, ICazaux &l Spaana (|2004l ) 
can account for a change of two orders of magnitude. How- 
ever, realistic theoretical estimations of the ^ factor would 
require simulations of galaxy evoluti on with spatial resolu - 
tion on the order of a few ~ 100 pc. (jCeverino et al.ll2010l ). 

Our approach to calculate the dust extinction is thus 



purely phenomenological. It does not assume any univer- 
sal extinction law for all galaxies and is not based on a 
dust production model. The extinction curve for each galaxy 
is different depending on its metallicity and gas contents. 
We can use the scaling between reddening, E{B — V), and 
extinction. Ay, to benchmark the impact of using a fixed 
universal extinction law. The C alzetti extinction cu rve has 
Rv = Av/E{B -V) ^ 4.05 (|Calzetti et al.ll2000l ). while 
for a Supernova (SN) extinction curv e it can be Ry = 7.8 
or Ry = 5.8 dHirashita et all 120051 ). Results at 2 ~ 5 
(|Forero-Romero et al. 2010t ) indicate that for the brightest 
best resolved galaxies {Muy < —21) Rv ~ 8.0 ± 0.5 is a fair 
approximation for the median values. This means that for 
a given amount of extinction, the expected reddening will 
be higher for a Calzetti extinction curve. In other words, 
one could match the UV magnitudes using the Calzetti ex- 
tinctio n law, as it was done for instance bv lDevriendt etahl 
(|2010l ), but conflicting results for the UV colours can be 
expected in that case. 

Our clumpy ISM model, as applied to the UV- 
continuum, fixes the optical depth of dust in the homoge- 
neous and clumpy phases of the IGM. The hydrogen optical 
depth in the homogeneous ISM is fixed by the HI column 
densities already calculated, while HI column densities in the 
clumpy phase can be bounded by the conditions expected 
in the molecular clouds of young star forming regions. With 
these constraints we proceed in the next section to quan- 
tify the expected extinction of the Lyman-a line based on 
a physical model that includes the resonance nature of the 
line. In order to be consistent with the approximation used 
for the continuum extinction we also fix the geometry of 
the gas and dust distribution to be that of a dusty slab with 
the radiation sources distributed homogeneously. In the next 
sections we will estimate the escape fraction of Lyman-a ra- 
diation in such configuration, and illustrate how this results 
are consistent with observational constraints of the escape 
fraction as a function of reddening. 



4 SLAB APPROXIMATION FOR LAES 

We approximate LAEs as homogeneous slabs of gas with 
dust and Lyman-a radiation sources homogeneously mixed. 
The motivation to explore such a model is twofold. First, 
we want to be consistent with the extinction approxima- 
tion already used for the UV continuum. Second, the ho- 
mogeneous distribution keeps an important feature seen in 
many simulations, including the best resolved galaxies in the 
MareNostrum Simulation used here, namely that the stars 
are not clustered around a single point with respect to the 
gas. 

Concerning the latter point, iLaursen et"all (|2009l ) stud- 
ied the escape fraction in high resolution simulations of in- 
dividual galaxies. The resolution study they performed in- 
dicates that converged values for the escape fraction require 
minimal smoothing lengths for the gas on the order of 160 

pc comoving, which is one order of magnitude smaller 
than the resolution of our MareNostrum High-z Galaxy For- 
mation Simulation. 

Therefore, it is still unavoidable to make use of the kind 
of subgrid models we propose here, in order to derive statis- 
tical results on the escape fraction. However, it is possible 
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Figure 1. Escape fraction for Lyman-« photons emitted inside 
different dusty gas configurations as a function of the product 
Tai^Tu)^^^ ■ Symbols represent the results obtained with CLARA. 
The empty triangles show the solution for a dusty sphere when the 
photons are emitted at the centre of the sphere, whereas stars rep- 
resent the case with Lyman-a sources homogeneously distributed 
in the sphere. Hexagons represent the results of the infinite dusty 
slab with sources homogeneously distributed. The solid line shows 
the analytical solution shown in Eq. ^ for the infinite dusty slab 
and sources located in the slab's centre = 0.525). The dashed 
and dotted line represent the same analytical expression with pa- 
rameters ^ = 0.625. A better fit for the homogeneous distributed 
sources (both for the sphere and the slab) is displayed by the 
dash-dotted lino, sec Eqs. JSj and ijOj. 



to improve the modeled physics by considering explicitly the 
effects of resonant scattering in the line absorption. 

In this Section we will employ CLARA to estimate the 
escape fraction in the slab configuration, assuming homoge- 
neously mixed sources. The source distribution constitutes 
the major difference between our work and similar Lyman-a 
Monte-Carlo radiative transfer studies (|Hansen fc Ohll2006l : 
IVerhamme et al1l2006l ). We will show that this assumption 
strongly affects the escape fraction results. We will then de- 
scribe how we estimate the relevant physical quantities in 
the simulation to obtain an escape fraction for each simu- 
lated galaxy. We show how this model agrees with observa- 
tional results on the escape fraction for galaxies at 2: ~ 2.2 
and z ~ 3. 



4.1 Radiative Transfer Results 



where tq is the Ifydrogen optical depth, Ta is the opti- 
cal depth of absorbing material (for albedo values of A, 
Ta — {1 ~ A)Td, where Td is the dust optical depth), and 
a is a measure of the temperature in the gas defined as 
a = Aul / {2Ai/d) , Apu ~ {vp/cjuo is the Doppler frequency 
width, and Vp — {2kT /mnY^'^ is -\/2 times the velocity dis- 
persion of the Hydrogen atom, T is the gas temperature, 
niH is the Hydrogen atom mass and Avl is the natural line 
width. The constant ^' = ^-/S/tt^''^^ and ^ is a free param- 
eter taking a value of ^ = 0.525 in the case of centrally 
located sources. 

For comparison, we show the results in the case of a ho- 
mogeneous dusty gas sphere with centrally located sources 
of Lyman-a radiation. This configuration is fully parame- 
terised by the optical depth of Hydrogen and dust as mea- 
sured from the centre of the sphere to its boundary. We 
simulate with CLARA a series of spheres with different values 
for these optical depths (see Table IXTj) . The empty triangles 
in Figure [T] represents the results for that configuration. The 
main point of this numerical experiment is to confirm that 
the same function that is used to describe the escape fraction 
of the infinite slab seems to be able to describe the escape 
fraction out of the sphere in the case of centrally located 
sources. The only difference is the value for the constant, in 
the spherical case ^ = 0.625. Nevertheless, the differences 
between both geometries for the gas density and the dust 
abundance, at a given temperature, are never larger than a 
factor of 4, for the range of parameters studied. 

We now turn to the case where the sources of Lyman-a 
radiation are homogeneously mixed inside the infinite slab. 
In this configuration we expect a higher escape fraction given 
the fact that the photons will be, on average, emitted closer 
to the escape surface. This will allow the photons to be af- 
fected by less collisions compared to the photons emitted at 
the centre of the slab. In Figure [T] the escape fraction for 
the homogeneous distribution of Lyman-a sources is shown 
by the hexagons. For small values of T = TaiaroY^^ the es- 
cape fraction seems to be well approximated by Eq.Q. For 
large values of T > 1 the escape fraction can be up to three 
orders of magnitude higher than in the central emission ge- 
ometry in the range of explored values. The scaling of the 
escape fraction with the variable T does not follow Eq.®. 
For large values of T > 1, the escape fraction shows a fall off 
almost proportional to 1/T. This can be understood in anal- 
ogy to the classical case of continuum light emitted inside 
a dusty slab. As the optical depth reaches high values, only 
the photons emitted close to the surface within a distance 
oc 1/r have a high escape probability. 

In analogy with the solution of continuum attenuation 
in a dusty slab, we find that for Lyman-a radiation, a solu- 
tion with the functional form 



The problem of an infinite homogeneous dusty gas slab has 
an analytic solution for the escape fraction provided that the 
sources are located in a thin plane. The dashed line in Fig- 
ure [1] corresponds to the theoretical expectation (described 
in Appendix 1X1 Eq. (|A19|l ) for the escape fraction out of a 
infinite slab as a function of the product {aroY^'^Ta, 



fesc — 



cosh (cV(a^o)i/^r„) 



(4) 



fesc — 



with 



1 - cxp(-P) 



e((aro) 



1/3^ ^3/4 



(5) 



(6) 



provides a reasonable description of the Monte-Carlo results 
as shown in Figure [T] (dot-dashed line) with e = 3.5 in the 
slab geometry and e — 1.0 in the spherical case. 

We conclude that the source distribution, relative to the 
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gas distribution, has a larger impact on the escape fraction 
than the geometrical distribution of the gas itself. The infi- 
nite slab exhibits very similar escape fraction compared to 
a spherical gas distribution, even the same scaling with re- 
spect to the gas properties holds. On the other hand, the ho- 
mogeneous source distribution presents a radically different 
scaling with the optical depth in the gas, allowing for high 
escape fractions at large values of {aroY^^Ta > 1, both for 
the slab and spherical geometries. Additionally, the shape of 
the outgoing spectra, corrected for the escape fraction are 
different for the homogeneous vs. central source distribution 
(see Appendix |X] for the case of the dusty sphere). 

4.2 Implementation in the MareNostrum 
Simulation 

In the previous subsection, we obtained results on the escape 
fraction for the case of a homogeneous dusty gas slab with 
radiation sources distributed homogeneously. We provided a 
fitting formula (Eq. (O) for the escape fraction as a function 
of the optical depths of gas and dust and the gas temper- 
ature. These results can be used by any galaxy evolution 
model that provides predictions for these quantities. 

In our case, we use this formula to to derive the values 
of the escape fraction for each galaxy in the MareNostrum 
High-z Universe simulation. In what follows, we explain how 
we proceed to estimate the physical quantities we need: tq, 
Ta and a. 

First, we recall from Section |3] that the dust model de- 
scribed above splits the extinction into two contributions: 

(i) the extinction by the homogeneous ISM on all stars in 
the galaxy, 

(ii) the extinction by the birth clouds around young stars. 

We now want to have estimates for the escape fraction 
associated with the homogeneous ISM, f'fc'', and with the 
birth clouds, f^sc ■ 

In order to estimate the ISM escape fraction, fif^^ we 
describe the ISM as a slab of constant density and temper- 
ature of 10*K, which fixes the value of a = 4.7 x 10""'. The 
hydrogen column density of this slab is given by Eq. ((2)1 . The 
average optical depth of neutral hydrogen is calculated from 
this column density and the hydrogen cross section at the 
center of the Lyman alpha line at a temperature of lO^'K. 
Based on the results of the previous subsection, we have thus 
determined all the parameters we need to calculate the es- 
cape fraction in the slab approximation: the optical depth of 
dust and Hydrogen (t|^*^, rlf"') and the gas temperature 
T. 

We now have to estimate the escape fraction associated 
with the birth clouds, f^^^. As all the Lyman-a emission 
comes from these regions, the full intrinsic Lyman-a lumi- 
nosity has to be corrected as well by this escape fraction. 
The dust optical depth associated with the clouds, r^*^, has 
already been calculated using Eq.Q. On the other hand, 
the optical depth of hydrogen in the birth clouds, r^"^, still 
has to be estimated. Based on observations of large molec- 
ular clouds, the neutral Hydrogen column density in the 
dense st regions is of the order of lO'^^ cm~^ jWannier et al.l 
1 19831 ). This is already ~ 3 orders of magnitude lower than 
the optical depth in our simulated galaxies at high redshift. 
Moreover, the warm regions have densities two orders of 



magnitudes lower l|Ferrierj|200ll 'l. A realistic estimate then 
puts the average neutral Hydrogen optical depth ~ 10^ 
times lower than the one estimated for the full galaxy. Un- 
der these conditions, we find within the cloud tq^ ~ r^'^ 
with ar^^ < 1 and t^'~^' > 1, making the extinction en- 
hancement by resonant scattering irrelevant. In this case, 
we can take the escape fraction as the continuum extinction 
at A = 1260 A for a spherical geometry. 

4.3 Comparison against observational constraints 

Observational estimates of the escape fraction at redshifts 
z > 5 (the lowest current redshift of the MareNostrum High- 
z Universe simulation) are not available. Nevertheless, at 
redshifts 2 ~ 3 the homogeneous ISM model is already con- 
sistent with the red dening scaling with galaxy luminosity 
(|Shaplev et al.ll200ll ). It is then reasonable to think that at 
all redshifts 2 < 3 the homogeneous ISM is still the appro- 
priate model regarding the extinction. For this reason, we 
decide to compare the escape fraction predicted by the ho- 
mogeneous ISM with the observational results at 2 ~ 2.2 
and 2 ~ 3. 

In the left panel of Figure [5] we present the redshift evo- 
lution of the escape fraction as a function of reddening at 
redshifts 2 ~ 5, 6, 7. We find that there is a weak redshift 
dependence. The biggest difference at each redshift is the 
presence of galaxies with larger values of colour excess. As 
the redshift evolution of the overall scaling is rather weak, 
we compare our theoretical results at 2 ~ 5, the last redshift 
available in the simulation to make this study, with the ob- 
servations at 2 ~ 2.2 and 2 ~ 3, the highest redshifts so far 
with observations of this kind. The physical time between 
the observational results and the simulated ones is ~ IGyr. 

In the right panel of Figure [2] we show the final re- 
sult of the comparison between the escape fractions from 
observations at 2 ~ 2.2 and 2 ~ 3 and our model ap- 
plied to the galaxies in the simulation at 2 ~ 5. The empty 
symbols represent the observational data, indeed showing a 
strong scaling: the escape fraction decreases with increas- 
ing colour excess. Our simple model of a dusty gas slab with 
homogeneously distributed sources reproduces w e ll the se ob- 
servational t r ends determined by iHaves et al.l (|2010l ) and 
iKornei et all l|2010l ). 

In order to facilitate the com parison of our re- 
sults with those from other works (|Kornei et al.l I2OI0I : 
iHaves et ahlbOlOl ). we also use the functional form fesc = 
CLyQlO"°*-^(«-^)*^!'° to fit the trend of the median values 
in the fesc-E{B—V) plane. This form takes into account that 
CLya 7^ 1 for very low extinction values {E{B — V) ~ 0) , as 
it is expected due to the resonance nature of the Lyman-a 
line. We obtain Clvc = 0.21 ± 0.05 and kLyc = 8.6 ± 1.0, 
which is s omewhat flatter than the value of kLya = 12.0 ob- 
tained by ICalzetti et al.l (|2000l ) for Lyman-o wavelengths. 



5 ANALYSIS AND RESULTS 

In this section we present our results on the escape fraction 
applying the slab model to the gas and dust contents cal- 
culated from the MareNostrum High-z Universe simulation. 
We first derive useful scalings of the escape fraction with the 
mass of the DM halo hosting the galaxy. We then apply the 
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Figure 2. Escape fraction as a function of the UV color excess as estimated from tiie liomogeneous ISM model. Tiie lines show different 
distributions obtained from the simulation, the thick line mar ks the median , and t he thin lines the first and third quartile. Stars represent 
observational estim ates of the escape fraction at ^ ~ 2.2 by iHaves et al.l 1 I2OICII ). Diamonds represent similar estimates at 2 ~ 3.0 by 
iKornei et al. 1 1 I2OIOI) In the left panel, the different lines show the evolution between 5 < 2 < 7. We observe how to lower redshift values, 
higher values of the color excess are present. The very weak redshift evolution justifies the comparison of the results of this model 
against observations shown in the right panel. Only the simulated galaxies above the resolution threshold were taken into account in 
the calculation. The slab model with homogeneous source distribution matches closely the observed scaling of the escape fraction with 
colour excess. 



estimated escape fraction to obtain the luminosity function 
of LAEs at high redshift, z > 5. Given that at these red- 
shifts the extinction model includes a component associated 
to the birth clouds of the stars, we will add this component 
for the comparison with the luminosity functions at these 
redshift s. 



5.1 Scaling with Halo Mass 

In the right panel of Figure [3] we show the total escape frac- 
tion Z™"^ — fife ' X f^J^ as a function of the host halo mass 
for the homogeneous slab model, where fi^J^' and /^J? were 
defined in sub- section 14.21 The overall scaling of fisJ^' and 
fesc with halo mass follows the relation parameterized by 
Eq.((S| where the product Ta{aToY^^ seems to be propor- 
tional to the halo mass. 

Due to this shape, the escape fraction has a large 
scatter for galaxies in haloes more massive than Mc = 
6 X lO^''/t~^M0. The escape fraction for haloes less mas- 
sive than this characteristic mass Mc seems to be always 
larger than 10%. 

We find that the median value of the escape fraction 
for the homogeneous ISM component, f'sJ^', between red- 
shifts 5 < 2 < 7 can be well approximated by the following 
expression: 



J esc 5 

where R depends on the halo mass, Mh, as follows: 



(7) 



R = 



Mh 
Mc 



5/3 



(8) 



and the critical mass is equal to Mc = 6 x 10^°/i"^Mq. 

Correspondingly, the median value of the escape frac- 
tion associated with the birth clouds of younger stellar pop- 
ulations, fUc, between redshifts 5 < 2; < 7 can be approxi- 
mated by the following expression: 



D 



where D is equal to 



20 Vm 



Mh 
Mc 



5/3 



(9) 



(10) 



and Mc has the same value as in the case of the homogeneous 
ISM. The results in the last two equations are dependent 
on the dust abundances derived by matching the observed 
UV luminos ity functions and the UV colours in the way 
presented in iForero- Romero et al.l (|201G| ). 

The observed scatter in /™"^ can be reproduced in 
a Monte Carlo fashion if one calculates individual escape 
fractions fif^' and fcs^ from halo masses in the simula- 
tion, using Eqs. ([7} and @ with different values of Mc and 
taking logj^f, i\fc//i^^MQ as a Gaussian variable with mean 
logio(6 X 10^") = 10.77 and variance of 0.3. 

In Figure 13] we compare our results against two other re- 
lated theoretical works at high redshift. We start by pointing 
out that the observed scaling in our results follows closely 
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Figure 3. Escape fraction as a function of darli matter iialo mass, for redshifts 2 ~ 5, 6 and 7. The tiiick lines mark tiie median, and the 
thin lines represent the first and third quartile. Triangles show different theoretical results. Left panel. The escape fraction includes the 
contribution from the homogeneous ISM jf l^r^^) only. The medi an values for f^gj^^ are reproduced by Eq. 10 The symbols correspond 
to the results for single halo simulations bv lLaursen et al Right panel. Contribution from the birth clouds around the younger 

population stars to the escape fraction (f^J^). The median values for are reproduced by Eq.(|9ll. The symbols represent the results 

obtained bv lPaval et al.l l|2010l ) from a cosmological simulation. 



the radiative transfer results shown in Figure [T] The corre- 
lation between halo mass and gaseous mass obtained in the 
simulation is translated here as a scaling between escape 
fraction and halo mass. 

In the left panel of Figure [3] we compare our results 
of the escape fraction from the ISM component against the 
results of high resolu tion simulations of individual galaxies 
(jLaursen et al.ll2009l ). The two models show a good agree- 
ment for high and low halo masses. There is, however, a dif- 
ference for haloes in the mass range ~ 3 x lO"/i"^M0. The 
mean of our escape fractions is higher by a factor of ~ 4, 
although the small sample of iLaursen et al.] (|2009h (seven 
galaxies) makes it difficult to estimate the statistical signif- 
icance of this comparison between the two models. 

In the right panel of Figure [3] we compare our results 
against a published model of the L yman-g escape fr action 
in a cosm ological cont e xt at 2 ~ 6 ( Daval et al.ll201ol '). The 
model of iDaval et"ai] (|201(]| ) presents a radically different 
prediction for the escape fraction in haloes more massive 
than Mth ~ lO"/i"^M0. While in our model the mean es- 
cape frac tion drops st e eply f rom ~ 0.1 down to ~ 10"'^, the 
results of iDaval et all (|201(]| ') show an increase of the escape 
fraction to values close to ~ 1 for the same range of mas- 
sive haloes. This disagre ement is very likely due to the fact 
that iDaval et all (|201(]| ) do not model explicitly the reso- 
nant nature of the line and approximate the ISM in such a 
way that the Lyman-a escape fraction is proportional to 
(1 — exp — {Td)) / Td where is the dust optical depth. This 
is a different approximation from the one that produces the 
results repr esented by Eq[5]and o ther similar radiative trans- 
fer studies (|Hansen fc Ohl l l2006f ) where there is an explicit 
dependence with the optical depth of neutral gas, thought to 



be abundant in the most massive and vigorous star forming 
galaxies at these redshifts. 



5.2 Intrinsic Emission and Luminosity Functions 

We now study the LAE LF predicted from our simulations 
taking into account the different contributions to the es- 
cape fraction. We will focus on the intrinsic emission associ- 
ated with star formation. The mechanism for star-triggered 
Lyman-a relates the amount of ionizing photons produced 
by young stars to the expected intrinsic Lyman-a luminos- 
ity. Here we assume that the number of Hydrogen ionizing 
photons per unit time is 1. 8 x 10^'' photons s~^ f or a star for- 
mation rate of 1 M0/yr (|Leitherer et al.lll999l ). This value 
assumes a Salpeter IMF, and that the star formation rate 
has been constant at least during the last 10 Myr. Assuming 
that 2/3 of these photons are converted to L yman-a pho- 
tons (case-B recombination. [Osterbroc^l 19891 ). the intrinsic 
Lyman-a luminosity as a function of the star formation rate 
is 



LLyc = 1.9 X 10* 



(SFR/Mq yr-^) erg s-\ 



(11) 



where the SFR is calculated from the total mass of stars 
produced in the last 200 Myr. We use this analytic expres- 
sion because the calculation of the total amount of ionizing 
photons directly from the information of the star particles 
in the simulation is not reliable given the short lifetime of 
the populations contributing to the ionizing flux and the 
limited resolution of our simulation. Only the star particles 
younger than 10 Myr will provide the bulk of the ionizing 
flux, these particles roughly correspond to 1% to 5% of the 
star particles in a given galaxy. The least resolved galaxy 
in our study has 200 star particles, which correspond to 2 
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to 10 particles to luUy sample the star formation history of 
the galaxy during the last 10 Myr. Only a handful of largest 
galaxies in the simulation (~ 100000 star particles, ~ 1000 
of them contributing to the ionizing flux) can give reliable 
results when the Lyman-a emission from the star particles 
is compared against Eq. (|lip . 

We find a tight correlation between star formation rate 
and halo mass for 5 < z <7. It can be approximated by 



SFR = 0.( 



Mh 



(lOio/i-iMo 



M: 



© yr 



(12) 



so that the final scaling between intrinsic Lyman-a emission 
and halo mass can be written as 



Lya 



1.29 X 10** 



lOiofe-iMo 



erg s 



(13) 



The observed Lyman-a luminosity can be calculated 
from our extinction model as: 



T° — ( fISM J old s , I rlSM pBC J- youngs 
J^Lya — {Jesc ^ ^Lya) "r VJesc ^ Jesc ^ Lya )' \^^) 

where the label old refers to the Lyman-a emission coming 
from stellar populations older than 25Myr and young labels 
the emission from stars younger than 25 Myr. Given that 
the intrinsic ionizing flux is negligible for stellar populations 
older than 25Myr, we can approximate the observed Lyman- 
a luminosity as: 



T ° 

J^Lya 



rISM rBC 

J esc J esc 



X I/Lya, (15) 

where the escape fractions fifj"' and /^J? are calculated as 
described in the previous subsections. 

Previous Monte-Carlo calculation s of Lyman-a radia - 
tive transfer in a multiphase medium (|Hansen fc OhI [20061 ) 
show that the dominant effect of the clumpy distribution of 
birth clouds on a photon propagating through the homoge- 
neous ISM is bouncing at the cloud's surface, justifying the 
approximation of not including further absorption by these 
clouds except at the Lyman-a source in the way we have 
just implemented. 

However, we note that the conditions for dust and gas 
abundance in our model are such that the calculations of 
iNeufeldl | 199lh for the escape fraction are not applicable 
here. The main assumption of that model (negligible absorp- 
tion and scattering in the homogeneous part of the ISM) are 
not met, meaning that the escape fraction in our model is 
not dominated by the effects of the clumpy component only. 

We have not taken into account the effects of scattering 
of Lyman-a photons by the IGM, after they escaped from 
the galaxy. This effect is very important before reionisation 
(2; > 6 in our simulation). After reionisation, it is expected 
that the blue part of the double peak spectrum would be 
absorbed as photons are cosmologically redshifted to reso- 
nance with neutral Hydrogen along the line of sight. To first 
order, we have taken this effect into account by dividing the 
intensity of the observed line strength I/^y^ by 2. 

In Figure |4] we show three different LAE LPs at three 
different redshifts 2 ~ 5, 6 and 7. The empty symbols repre- 
sent the observed LPs. The observational results at 2; ~ 4.5 
come from 50 spectroscopically confirmed LAEs from the 
Large Area Lyman-a survey (LALA) which covered a field 



of ~ 0.7 deg^ correspo nding to a comovin g survey vol- 
ume of 7.4 X lO' ^ Mpc^ (iDawson et alll2007l '). At z ~ 4.86 
the results from IShiova et al.l ( 20091 ) are based on observa- 
tions of the Cosmic Evolution Survey (COSMOS) field of 
1.83 deg'^ giving an effective volu me of 1.1 x 1 0^ M pc'^. 
At z ~ 6 the data are taken from lOuchi et al.l (|2008l ). m 
this case the LP was estimated from LAEs from 1 deg^ 
Subaru/XMM-Newton Deep Survey (SXDS), which probed 
a comoving volume of ~ 10^ Mpc^ . We include as well data 
at 2 ~ 5.7 in a field of 1 . 95 de g^ covered by the COSMOS 
survey iMuravama et all (|2007l ') and recent estimations at 
4.5 < z < 6.6 from the Vimos-VLT Deep Survey (VVDS) 
jCassata et al.|[201ll ) . The LP at z ~ 7 was constructed from 
a sample of 17 LAEs confirmed spectroscopically and 58 pho- 
tometr ic LAEs probing a com oving volume of ~ 2.17 x 10^ 
Mpc^ (|Kashikawa et al.ll2006h . In our simulation we probe 
a comoving volume that is of the same order of magnitude 
of these surveys. 

In Figure |4] we represent with black squares the LP 
estimated from the intrinsic luminosities calculated from 
the star formation rates. The spatial abundance is already 
corrected from the different halo abundance between the 
WMAP5 and the WMAPl cosmologies. Compared with the 
observational estimates at all redshifts, the normalisation is 
at least one order of magnitude higher. 

We can now apply the expected correction by the es- 
timated escape fraction of each galaxy. In our model there 
is a strong correlation between the mass of the halo hosting 
the galaxy and its associated escape fraction. In a simpli- 
fied manner, galaxies in haloes with masses smaller than 
6.0 x 10^° /i~^MQremain practically unaffected. This is seen 
in the LPs in Figure |4] where the hexagons represent the 
luminosity function corrected by the escape fraction. The 
bright end of the LP is modified but the overall normalisa- 
tion is kept one order of magnitude higher than observations. 

This result is completely analogous to the situation 
found in the continuum U V using the same simulation 
(|Forero- Romero et al.|[201ol ). The strong scaling of the red- 
dening with galaxy mass causes only the bright end of the 
UV LP to be effectively modified when extinction is taken 
into account. In the context of our model, applying a con- 
stant reddening value to all galaxies can only be physically 
justified if there is an additional extinction term from the 
youngest stars, as described in Section [S] 

If we add now the expected correction of the birth 
clouds around stellar populations younger than 25 Myr, we 
find that the faint end of the LAE LP can be further modi- 
fied. However, at 2 ~ 5 and 2 ~ 6 there remains an excess at 
the faint end of the LP, corresponding to haloes with masses 
1.0 X lO"fe"^M0 < Mh < 4.0 X lO^°/i"^M0. This result 
suggests that either the escape fraction for these galaxies or 
their star formation rate is too high. Nevertheless, the bright 
end of the simulated LPs, both at z ~ 5 and 2 ~ 6, shows a 
very close agreement with the observed one. The results are 
within Poissonian uncertainty. We have also checked that 
if we apply the scalings for the intrinsic emission and es- 
cape fractions with halo mass, Eqs. (f71[51[T5|). to a pure DM 
only simulation with a cubic volume of 250 /i^^Mpc on a 
sid e and 2048^ particl es (the Bolshoi Simulation presented 
bv lKlvpin et al.l (|2010l )') split into smaller sub-boxes, we can 
reproduce the results we have just discussed. Specifically, the 
scatter at the bright end due to cosmic variance is consistent 
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Figure 4. Luminosity Functions at 2: ~ 5, 6 and 7. The observational constraints (empty symbols) are compared to different results 
from tfie simulation (filled symbols): Intrinsic Lyman-a emission considering the total star formation rate in Eg.l lllll (black squares), 
Lyman-o emission corrected only by the extinction of the homogeneous ISM based on the spherical model studied with the Monte Carlo 
simulations using Eq. (O (blue triangles), Lyman-« emission corrected only by the extinction of the birth clouds using Eq. l[9ll (green 
hexagons). The final luminosity function corrected by the homogeneous ISM and the birth clouds is shown by the red diamonds. 



with the observations. Nevertheless, the same scatter does 
not help to explain the lower abundance of our numerical 
LAEs at the faint end. 

The normalization of our LF functions at redshift 2 ~ 7 
is still higher than observed. In principle, it could be possible 
to account for this difference by properly modeling the IGM 
absorption, which at this epoch, it is not completely ionized 
and will add a dimming effect to the LAEs. Considering a 
dimming factor that depends on luminosity only (without 
any scatter), we would require that 25% (40%) of Lya pho- 
tons be transmitted through the IGM in order to match the 
bright (faint) end of the observed luminosity function. Using 
again the DM only simulation, together with the scalings 
we obtain for the star formation and the escape fraction, 
we find that cosmic variance can account for a 0.5 dex (0.3 
dex) variations at the bright (faint) end of the luminosity 
function at z ~ 7. In conclusion, applying a constant IGM 
transmission of T = 0.3 — 0.4 we can still reproduce the lu- 
minosity function for the brightest three bins in Figure [4] A 
m ore realistic treatm ent of the effects of IGM (in the spirit 
of I Zheng et all (|2010l )) actually yields a large scatter for the 
transmission at a given Lynian-Q luminosity. Nevertheless, 
a detailed modeling of the IGM effect is far beyond the scope 
of this paper. 



6 CONCLUSIONS 

In this paper we model the escape fraction of Lyman-a pho- 
tons in the approximation of a dusty gas slab with Lyman-a 
sources homogeneously mixed. The escape fractions for this 
configuration and different dust and gas contents have been 
calculated using CLARA, our new Monte-Carlo code described 
in detail in AppendixlX] These results can be applied in any 
model that predicts the optical depth of gas and dust in 
galaxies. 

We selected the slab geometrical configuration in or- 
der to be consistent with the assumptions that led us to fix 



the dust abundances from the MareNostrum High-z Uni- 
verse simulation, as constrained by high reds hift UV at 
2: > 5 observations jForero-Romero et al.|[201ol ). Our pro- 
posed dust model describes the contributions from an homo- 
geneous ISM (slab geometry) and a clumpy phase (spherical 
geometry) associated to stellar populations younger than 
25Myr. 

We estimate the scaling of the Lyman-a escape frac- 
tion with the expected reddening for the slab component 
and dust abundances in the MareNostrum High-z Universe. 
The scaling shows a weak redshift dependence between 
5 < 2 < 7, furthermore there is a very good agreement with 
the observational estimation of the escape fraction and red- 
dening for galaxies at z ~ 2.2 and 2 ~ 3 assuming that only 
the homogeneous ISM component is dominant at these red- 
shifts. Including both contributions from the homogeneous 
ISM and the clumpy phase, we have calibrated the escape 
fraction as a function of the host dark matter halo based on 
the results of the MaraNostrum High-z Universe. 

As an application of these results, we construct the in- 
trinsic LAEs LF as estimated from the star formation rates. 
We find that the normalisation is one order of magnitude 
higher than observational estimates. Correcting the intrinsic 
LAE luminosities by the estimated escape fraction (homoge- 
neous ISM and clumpy phase included) brings the simulation 
into agreement with the observations at z ~ 5 and 2 ~ 6. 
The match at the bright end is acceptable within the Pois- 
sonian and cosmic variance errors. The mismatch at 2 ~ 7 
can be explained because a proper modeling of these epochs 
has to account for the yet incomplete reionisation process 
and the infiuence of the neutral parts of the IGM. 

Nevertheless, our results are in conflict with the ob- 
servational estimates at the faint end. There seems to be 
an excessive production of LAEs with intrinsic luminosities 
Ltd ~ 6.0 X 10'''^ erg s~^ at 2 ~ 5 and 2 ~ 6. This excess 
wa s spotted, although weakly, in the results of the UV LFs 
of (|Forero- Romero et all [20101 ). Furthermore, a fine tuning 
of the extinction model for galaxies in this mass range could 
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provide a better match to the observations. In that case, the 
galaxies at the faint end should be more dusty, which would 
make them redder, and thus breaking the broad agreement 
for the UV colours that we have already obtained. This is 
not a satisfactory approach to explain for these differences. 
Based on these considerations, we think that a probable ori- 
gin of this discrepancy is the high rate of star formation in 
galaxies situated in these haloes. A possible physical expla- 
nation is that supernova feedback modulates more effectively 
the star formation in haloes of masses < lO^^/i^^M©, pro- 
viding a mechanism to shape the faint end of the luminosity 
function. However, it is also possible that the star formation 
rate is overestimated in all haloes at these redshifts. This 
would be translated into a trade-off with less dust extinc- 
tion to explain the UV luminosity function. Regardless of 
what is the correct explanation of this enigma, all possible 
solutions seem to challenge our current understanding of the 
rate at which gas is converted into stars at high redsliift. 

It is encouraging that the results for the brightest galax- 
ies, the best resolved ones, both in observations and simu- 
lations, are consistent with observations in the UV (mag- 
nitudes and colours) as well with the Lyman-a line. A full 
picture of these massive high redsliift galaxies will be com- 
pleted by observations of the rest frame IR to be performed 
by ALMA. We will address in a upcoming work the predic- 
tions of our model in terms of ALMA observations, focus- 
ing on the most massive galaxies, constructing a complete 
panchromatic perspective of high redshift galaxies in the 
MareNostrum High-z Universe. 
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APPENDIX A: LYMAN-a RADIATIVE 
TRANSFER 

The basic principle of CLARA is to follow the individual scat- 
terings of a Lyman-a photon as it travels through a distribu- 
tion of gaseous hydrogen. Each scattering, which is in fact an 
absorption and re-emission, does not modify the frequency 
of the photon in the rest-frame of the hydrogen atom. But 
due to the peculiar velocities of the atom in the new di- 
rection of propagation of the photon, the new frequency in 
a laboratory rest-frame is different from the incoming fre- 
quency. Thus the photon performs a random walk not only 
in space but also in frequency. 

The relevant properties of the gaseous hydrogen can be 
fully described by its density, temperature and bulk veloc- 
ity. This is sufHcient to describe the emergent spectra of 
a source of Lyman-o? photons embedded in gaseous hydro- 
gen. It is important to observe that none of the emitted 
photons is completely lost by absorption as it is immedi- 
ately re-emitted. The original spectrum morphology of the 
Lyman-a source is modified, but the only way to lose the 
energy input by a Lyman-Q source is through dust absorp- 
tion. A simple description of the dust abundance in the gas 
must then be included to calculate its effect on the outgoing 
properties of the traveling Lyman-a photons. 

In the next subsections we give a detailed account on 
the basic underlying physics of the qualitative description 
we have just given. Once the physical fundamentals are de- 
scribed, we describe how these are implemented in the code. 
The last subsection is devoted to show the results available 
analytical test cases we applied the code to. 
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Al Physical Principles 

A 1.1 Scattering 

The scattering cross section of a Lyman-a photon is a func- 
tion of the photon frequency. In the rest-frame of the atom 
it is equal to 

where /12 = 0.4162 is the Lyman-a oscillator strength, vo = 
2.466 X lO^^Hz is the line centre frequency, Aul = 4.03 x 
lO^^uo = 9.936 X 10^ Hz is the natural line width, and the 
others symbols conserve their usual meaning. 

In the case of a Maxwellian distribution of atom veloc- 
ities, after convolving the individual cross section with the 
atom velocity distribution, we can write down the average 
cross section as 



W{e) ocl + - cos^ e (A5) 

where R/Q = for the 2P1/2 — >■ lS'1/2 transition and R/Q — 
3/7 for the 2P3/2 — >■ IS1/2 case. 

The spin multiplicity of each sate is 2J -f 1, meaning 
that the probability of being excited to the 2P2/3 state is 
twice that of the 2P1/2 state. These r esults are valid only 
in the core of the line profile. IStenfiol (|l976l ) found that in 
the wings quantum mechanical interference between the two 
lines act in such a way as to give a scattering resembling a 
classical oscillator. In that case the phase function takes the 
form 

w{e)xi + cos^ e (A6) 

The travelling distance, of a Lyman-a photon of fre- 
quency X can be expressed as 



a{x) = ax ^ fi2- 



where -ff (a, 1) is the Voigt function. 



_ff (a, x) 



H{a,x), 



(A2) 



(a- — + a- 



(A3) 



AfD = {vp/cjfo is the Doppler frequency width, Vp — 
{2kT /muY^^ is y/2 times the velocity dispersion of the Hy- 
drogen atom, T is the gas temperature, mn is the Hydrogen 
atom mass, a — Az^l/(2Ai^d) is the relative line width and 
x = {vi — vo)/Ai'D is a re-parameterisation of the photon 
frequency respect to the line centre normalised by the tem- 
perature dependent Doppler frequency width of the gas. 

The scattering is coherent in the rest-frame of the pho- 
ton, but to an external observer, any motion of the atom will 
add a Doppler shift to the photon. Measuring the velocity 
of the atom, Va, in units of thermal velocity u = -Va/vp, the 
frequency in the reference frame of the atom is 



x' = X — u ■ Ui, (A4) 

where rii is a unit vector in the incoming direction of the 
photon. 

In the general case, the scattering of the Lyman-a atom 
is not isotropic. For symmetry reasons the scattering is 
isotropic in the azimuthal direction, respect to the outgoing 
scattering direction. The distribution of the scatter direc- 
tions depends only on the angle, 0, between the incoming 
and outgoing direction of the photons, and rio, respec- 
tively. 

The information on the outgoing angle is encoded in 
the phase function, W{6). In general the angular momenta of 
the initial, intermediate and final states are involved in the 
calculation of W{6). In the case of resonant scatter the initial 
and final state are the same. The intermediate state corre- 
sponds to the excited state. Following the notation n/ jThere 
are two possible excited states, 2Piy2 and 2P3/2- 

Specifically, in the dipole approximation the phase func- 
tion can be written as: 



(A7) 



where nn is the neutral hydrogen number density, and it 
has been assumed that along the path I the temperature 
and bulk velocity field of the gas are constant, to ensure 
that the photon frequency can be represented by the same 
value of X along its trajectory. 

In what follows, we will always characterise an homoge- 
neous and static medium using the optical depth ro at the 
line center. 



A1.2 Dust Absorption 

In the case of dust interaction the photon can either be 
scattered or absorbed. The optical depth of dust, ad, can be 
generally expressed as the sum of an absorption cross section 
(Ja, and a scattering cross section as- 



ad ~ aa + as 



(A8) 



The determination of these two cross sections can be 
achieved using two different approaches. We name the first 
approach the ah initio approach. The ab initio approach 
seeks to determine the values of the dust cross sections from 
individual studies of dust grain properties and its interac- 
tion w ith photons. This is the approach used in the stud- 
ies by IVerhamme et all (|2006l ) . The second approach is a 
phenomenological one and it defines the dust properties in 
relation to the gas properties in the galaxy. The dust cross 
section properties are then derived from observations. In the 
interest of keeping our model simple to operate and with a 
good match to the level of detail required for the approxi- 
mation, we proceed with an ab initio approach. 

We express then the absorption and scattering cross 
sections as 



— nd 



(A9) 



where d is represents an average dust radius and Qa,s is 
an absorption/scattering efficiency. The dust albedo can be 
then expressed as 
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(AlO) 



At UV wavelengths the emission and absorption pro- 
cesses are equaUy hkely, with Qa ~ Qs ~ 1, making the 
dust albedo around ~ 1/2. We now express the dust of op- 
tical depth, Td in an analogous way to the neutral hydrogen 
optical depth Ea. (|A7P for a parcel of dust of linear dimen- 
sions I 



A2.2 Photon Displacement and Interaction 

For each photon of frequency x and direction of propagation 
Hi, we determine the distance freely travelled until the next 
interaction. This optical depth is determined by sampling 
the probability distribution Pij) = exp(— r) by setting Ts = 
— \n{TZ), where 7?. is a represents a random number between 
< 7?. < 1 following an uniform distribution. 

This optical depth interaction fixes the travel distance 
Is by: 



Td = O-dUdl, 



(All) 



Is 



(A12) 



where Ud represents the number density of dust particles 
and it has been assumed again that the dust cross section 
and dust number density are constant on the scale of I. 

A2 The Radiative Transfer Code 

The code implements a Monte Carlo approach to the ra- 
diative transfer. CLARA follows the successive scattering of 
individual photons as they travel through the gas distribu- 
tion, changing at each scatter the direction of propagation 
and frequency of the photon. We describe now in detail the 
technical implementation. 

A2.1 Initial Conditions 

The problem to solve defines the physical characteristics of 
the gas distribution and the initial conditions of the Lyman- 
Q emitted photons. The gas is described by the following 
characteristics: 

• The geometry of the gas distribution. In this section we 
will present results for the following configurations: infinite 
slab and sphere 

• The size of the gas distribution. This is parameterised 
by the hydrogen optical depth, to at line centre a; = 0. In 
all the geometries we measure the optical depth from the 
centre of the configuration to its nearest border. 

• The temperature of the gas distribution, T. This is set 
to a constant for all the gas distributions explored in this 
paper. 

• The gas bulk velocity field, Vb(r). This is in general 
dependent on the position. The only bulk velocity field ex- 
plored in this work corresponds to a Hubble like fiow in the 
spherical geometries. 

• The dust optical depth, Td- 

• The dust albedo, A. 

The photons are described by the following properties: 

• The spatial distribution with respect to the gas dis- 
tribution. There are two possibilities. All the photons are 
emitted from the centre of the gas distribution or they are 
homogeneously distributed throughout the gas volume. 

• The initial direction of propagation. We assume that 
the emission is isotropic (in the local comoving frame). 

• The initial frequency, Xi. We consider that all photons 
are injected in the centre of the line Xi — 0. 

The number of photons used in order to reach con- 
vergence are Nph ~ 5 x 10''. In some cases convergence is 
reached with Nph ~ 1 x 10'' photons. 



The photon travels a distance h in the direction bi, 
at this point we decide if the photon interacts either with 
an hydrogen atom or with a dust grain. The probability of 
being scattered by a Hydrogen atom is 



Ph = 



o-xnH 



o-xnH + OdUd 



(A13) 



Another random number TZ is generated and compared 
to Ph- \i TZ < Ph the photon interacts with Hydrogen, 
otherwise it interacts with the dust grain. In the case of 
interaction with dust a new random number is generated 
and compared with the albedo A. If the random number is 
less than the albedo, the photon is scattered coherently in 
a random direction, otherwise the photon is absorbed and 
not considered for further scatterings during the rest of the 
simulation. 

In the case of photon interaction, the situation is more 
involved. As already described by Eq. (|A4|) . the new pho- 
ton frequency in the observer frame depends on the specific 
velocity of the atom, u. It is therefore necessary to draw a 
velocity for the hydrogen atom. The preferred velocity for 
the atom cannot be drawn from a isotropic Gaussian distri- 
bution, there is an implicit spatial anisotropy in this case. 
The opacity of the gas is frequency dependent: the preferred 
atom velocity in the parallel direction to the photon prop- 
agation is different from the velocity in the perpendicular 
direction. 

The atom velocity is thus determined in two steps. In 
the first step the two perpendicular velocities are selected 
from a Gaussian distribution. In the second step the paral- 
lel component of the atom velocity is calculated. The par- 
allel component is calculated from a distribution calculated 
as the convolution of a Gaussian (representing the intrinsic 
velocity) convolved with the Lorentzian probability of the 
atom being able to scatter the photon: 



2 a 

/(ujl) oc Q{u\\) X C{x — M||) oc e "II X 



■K {x — ui^\)'^ -\- 



The resulting normalised probability is: 

2 

n ^ W 



TTH{a,x) {x — u\\)^ + 



(A14) 



The distribution represented by Eq. (|A14[I is not analyt- 
ically integrable, therefore we generate the parallel velocities 
M|| by means of the rejection method. 

We stress that we do not use any of the commonly used 
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speeding mechanisms. These speeding techniques are moti- 
vated by the physical fact that scattering of atoms in the 
core frequencies are irrelevant from the point of view of the 
change in frequency and position that they represent. Given 
the uncertainty of applying these techniques (despite of care- 
ful code calibration) in the presence of dust, we decide to 
use the simplest rejection method. 

The new emission direction of propagation (in the atom 
rest-frame), n/ is determined as already described in detail 
in subsection I All In the observers frame, the final frequency 



X f ^ Xi — Ui ■ u + Uf ■ u + g{l — ni ■ n/), (A15) 
where the factor g accounts for the recoil effect. 



A3 Analytical Tests 

The basic tests on CLARA are based on the available analyt- 
ical solutions for a thick, plane-parallel, isothermal infinite 
gas slab of uniform density. The solution of the infinite slab 
provides the analytical expressions for: the distribution of 
outgoing frequencies, the average number of scatterings and 
the position of the peaks of the outgoing spectrum. In the 
case of a dusty infinite slab one we also have an analytical 
expre ssion for the escape fraction (|HarrinErtonlll973l : lNeufeldl 

S). 

The expression for the analytic emergent spectrum as a 
function of frequency shift for a mid-plane source reads: 



J{x) = 



24 aro cosh(V7r4/54|2;3 - x^\/aTo) ' 



(A16) 



where to is the optical depth measured from the mid-plane 
to the slice boundary and Xi is the frequency of the photons 
injected in the slice. For Xi = the emergent spectrum has 

a symmetric two-peaked shape. 

Setting dJ/dx = (|Harringtonll 19731) . it can be shown 
that the position of the maxima is: 



x^ = ±1.066(aro)'''^ 
while the average number of scatterings is: 

(N,) = 1.612ro. 



(A17) 



(A18) 



If dust is added to the gas distribution, it is then pos- 
sible from analytical considerations to estimate the scalings 
of the escape fraction from the slab physical properties. The 
expression for the escape fraction reads 



(A19) 



cosh 



(e\/(^^ 



where ^' = ^VS/tv'^^'^^ and ^ = 0.565 is a parameter found 
from a fit. 

In the following we compare the analytical results de- 
scribed by Eqs. (|A16|) to HA19P with the results from CLARA . 
In Figure \KT\ we show the outgoing spectra for the infinite 
slab with four different optical depths tq = 10^, 10®, 10^ 
and 10*, and a constant temperature of T = lOK. The broad 



shape of the double peak is reproduced. The agreement with 
the expected analytical solution becomes better as the hy- 
drogen optical depth increases, as expected. The simulation 
and the analytical result are expected to hold for optical 
depths and temperatures such as aro > 10'^. 

In the middle panel of Figure IXTl we have the positions 
of the positive and negative peaks in the outgoing spectra, 
as a function of the product aro. These results were ob- 
tained from 7 different runs with their characteristics listed 
in Table lAll In the right panel we compare the number of 
scatterings with the theoretical prediction, which is only de- 
pendent on the hydrogen optical depth tq. In both cases we 
find optimal agreement with the theoretical expectations. 
We stress again that even though we reproduce the right 
number of scatterings, we have not implemented any accel- 
eration scheme based on the skip scattering scheme. 

In Figure IA2I we show the results concerning the es- 
cape fraction for the dusty infinite slab. We compare our 
results to the expected result from Eq. (|A19|) . In Table \KT\ 
we list the parameters for the dusty slab run, as well as the 
values of the escape fraction we find. The theoretical expec- 
tations for the escape fraction take the interaction with dust 
as an absorption. In the case of our runs, as we consider an 
albedo of ^ = 0.5, the optical depth of dust grains has to 
be replaced by an effective value of an absorbing material 
Ta = {1 — AjTd = 0.5rd, where r^ is defined by Eq. (lAllf) . 
Once again, we find a reasonable agreement with the theo- 
retical expectations. 

At this point we are confident that CLARA provides the 
expected outgoing spectra and number of scatterings. The 
dust implementation has been tested against the available 
theoretical constraints, yielding satisfactory results. 

The case of sphe rical symmetry offe rs analytical solu- 
tions in specific cases (|Diikstra et al. [20061 ). Besides, the case 
of a expanding (contracting) sphere has been extensively 
studied in the literature with different Monte Carlo codes. 
The spherical symmetry situations gives us the possibility to 
perform further tests on CLARA. We will deal with the case 
of a static, isothermal and homogeneous sphere of gas first. 
Sub sequently we w il l add a Hubble-like flow to the sphere. 

iDiikstra etlo] (|2006l 'l computed the analytic solution 
for the emergent spectrum of a point-like Lyman-Q source 
surrounded by an homogeneous, static gas distribution at a 
constant temperature. The solution can be expressed with 
the same functional form as the H arrington-Neufeld solution 
(|Harringtonlll973l : lNeufeldlll99ll ). 

We setup different configurations in CLARA for a series of 
optical depth and temperatures. Figurc [A3l shows CLARA out- 
puts compared with the corresponding analytical solutions. 
Once more, in the limit of large values for arc > 10^ we 
recover the behaviour exp e cted by the analytical solution 
provided bv IDiikstra et al.l (|2006l ). 

We add to the homogeneous sphere a radial Hubble like 
flow as a function to the distance to the centre of the sphere 



Vr{r) = Vmax-f:, 



(A20) 



where Vmax is a constant which can be selected to be positive 
or negative. 

The outgoing spectra of this type of configuration has 
not been worked out from an analytical point of view. How- 



16 Forero- Romero et al. 
15 



10 



1 1 1 1 1 1 1 1 

TO =10= 


1 


1 1 1 1 1 1 1 1 


__- T0=10« 






TO =10^ 






-_ —- T0=108 




- 


- li 








If 




-50 





50 




X 





2.0 



1.5 - 



0.5 



T — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — r 

■ Harrington (1973) 

■ Negative Xm 

Positive Xm 



J I I I I I I I I I I I I I I I L- 

2 3 4 

logio aTo 



8 - 



I I I I I I I I I I I I I 



T-r-r 



Harrington (1973) 
This Worl< 



^ I I I I I I I I I I ■ ■ ■ ■ I I I 

5 6 7 8 

logio To 



Figure Al. Analytical tests related to the infinite homogeneous gas slab. The left panel show the outgoing spectra for different hydrogen 
optical depths, compared to the expected analytical solutions. The middle panel compares the positions of the maxima in the outgoing 
spectra to the expected analytical estimates for a series of 5 different runs. The panel on the right compares the average number of 
scatterings the Lyman-« photons go through before escaping the slab. CLARA passes successfully these three tests. 
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Table Al. Summary of the physical characteristics for all runs, both with infinite slab and spherical geometry of a dusty gas distribution. 
Vmax refers to the parameterisation of the radial velocity profile imposed on the sphere following a Hubble-like law. The escape fraction 
in the case of centrally distributed sources in a sphere is denoted by f '-^ , , while . corresponds to the escape fraction when the 

Lyman-Q sources are homogeneously distributed. The escape fractions for the infinite slab are denoted by ^^^^ and slab' 



ever, it is part of the most studied physical setups by Monte 
Carlo codes similar to ours. The quantitative behaviour of 
our results follows the behaviour of previously published re- 
sults. In the expanding case, there is suppression of the blue 
peak accompanied by an enhancement in the red wing of 
the emergent spectrum. The case of the collapsing sphere 
displays the enhancement in the blue peak and suppression 
in the red peak. 

A4 The Homogeneous Source Distribution 

The model of the LAEs suggested in the paper is based 
on the approximation of an homogeneous dusty gas sphere, 
with the photon sources distributed homogeneously inside 



the sphere. In the body of the paper we explored in detail 
the consequences of this model on the escape fraction. At 
the same time we have studied the effect on the outgoing 
spectra. 

Figure IA5I presents the outgoing spectra already cor- 
rected by the escape fraction. The dashed line shows the 
result for the sphere with sources located at the center, the 
solid line is the expected result for the case of the dustless 
sphere. The match between the two curves shows that ef- 
fectively the shape of the outgoing spectra is preserved in 
the presence of dust. This result is expected because all the 
photons are emitted at the same location and experience the 
same optical depth. 

In case of homogeneously distributed sources, the spec- 
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Figure A2. Escape fraction for Lyman-o photons emitted from 
the centre of a infinite slab of gas. The hexagons represent the 
results of CLARA, the continuous Unc, the analytical solution of 
Eo. lTSlQt . 



10 



1 — I — I — I — I — I — I — I — I — I — I — I — r 

Dijkstra et al.(2006) 

1 = ^0l 

x = 10^ 

x = 10^ 




T 



T 



T 



— Dijkstra et al.(2006) 
. -- Static 

■ ■■ Expanding 

— Contracting 




20 40 



Figure A4. Emergent line profile in the coUapsing/outfiowing 
sphere tests. The collapsing/outflowing solutions arc symmetric 
under parity transformations x —x. 
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Figure A3. Emergent line profile for the static and dustless 
sphere test. The photons are emitted at the centre of a sphere 
with d ifferent hydrog e n opt ical depths. The analytical solution 



with d iiierent hydrog e n opt i 
follows [Dii"kstra et al.l ||2006|) 



tral shape is different. The peaks are located closer to the 
line centre due to the fact that now, on average, the photons 
experience a smaller optical depth. At the same time, each 
wing presents now an asymmetrical shape with respect to 
these maxima. Furthermore, there is now some amount of 
photons very close to the center of the line, i.e. emitted close 
to the surface, that escaped without any scatter. 

We have presented CLARA, a radiative transfer code for 
Lyman-a radiation. The code propagates Lyman-a photons 
through arbitrary distributions of hydrogen density, temper- 
ature, velocity structures and dust distributions. The code 
has successfully passed all the available tests for codes of its 
kind, and has allowed us to study simplified geometries in 
the present paper. 



Figure A5. Emergent line profile in the case of a static dusty 
sphere. The parameters for the run where T = lO^K, tq = lO'^ 
and Tjj = 0.4. The profiles are corrected for the escape fraction. 
The different line profiles correspond to the case of centrally dis- 
tributed sources versus homogeneously distributed. The case of 
centrally distributed sources i s compatible with the analytic so- 
lution of iDiikstra et ah The case of homogeneously dis- 
tributed sources presents a different line profile. The maxima are 
closer to the line center, and each peak is asymmetric respect to 
the maximum. 



